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Bifurcations leading to complex dynamical behaviour of non-linear sys- 
tems are often encountered when the characteristics of feedback circuits in 
the system are varied. In systems with many unknown or varying parame- 
ters, it is an interesting, but difficult problem to find parameter values for 
which specific bifurcations occur. In this paper, we develop a loop breaking 
approach to evaluate the influence of parameter values on feedback circuit 
characteristics. This approach allows a theoretical classification of feedback 
circuit characteristics related to possible bifurcations in the system. Based 
on the theoretical results, a numerical algorithm for bifurcation search in a 
possibly high-dimensional parameter space is developed. The application 
of the proposed algorithm is illustrated by searching for a Hopf bifurcation 
in a model of the mitogen activated protein kinase (MAPK) cascade, which 
is a classical example for biochemical signal transduction. 

1 Introduction 

A frequent challenge in the analysis of non-linear dynamical systems is to find pa- 
rameter values for which the system undergoes changes in its dynamical behaviour. 
Such changes are directly related to the emergence of complex dynamical behaviour. 
Standard cases of complex dynamical behaviour are multistability, i.e. the existence of 
several stable steady states, limit cycle oscillations, and non-periodic oscillations. 
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Feedback circuits are the major structural feature in the emergence of complex 
dynamical behaviour. In particular, it can be shown that a positive feedback circuit in 
the system is required for multistationarity , whereas a negative circuit is typically 
required for limit cycle oscillations [1^. This importance of feedback circuits makes 
control theory a natural tool for the analysis of complex dynamical behaviour. 

Yet, the main properties of a system's qualitative dynamical behaviour are the loca- 
tion and stability of equilibrium points. Knowledge of these is often also useful when 
analysing complex dynamical behaviour. It is well known from dynamical systems the- 
ory that two stable equilibrium points are separated by an invariant repellor, which 
contains an unstable equilibrium point in most cases. Similarly, stable limit cycle os- 
cillations usually coexist with an unstable equilibrium point. Also transient behaviour 
is often governed by the attraction to and repulsion from equilibrium points. Thus a 
convenient first step when studying the qualitative behaviour of a dynamical system 
is to look at stability properties of equilibrium points. 

A classical tool for analysing the influence of parameter values on the location and 
stability of equilibrium points is bifurcation analysis. Bifurcation analysis is done rou- 
tinely with numerical continuation methods for one adjustable bifurcation parameter 
[l^ . Methods for numerical bifurcation analysis in several parameters are now being 
developed (sl. [27|. but due to practical considerations, they remain limited to two or 
three adjustable bifurcation parameters. 

The challenge to find parameter values for bifurcations is of particular relevance in 
the area of biological systems. The main reasons for this are that biological function 
is often based on complex dynamical behaviour, and that parameters can vary within 
a large range due to environmental or internal conditions. 

There are many examples where complex dynamical behaviour of a non-linear bi- 
ological system can directly be related to biological function. Some examples from 
the specific area of biochemical signal transduction within living cells are given by 
bistability in the mitogen activated protein kinase (MAPK) pathway to induce devel- 
opmental processes , rapid activation of caspases upon an over-threshold stimulus 
in programmed cell death y, and sustained oscillations in circadian clocks 13l. |22|. 

Systems for biochemical signal transduction are usually modelled with non-linear 
ordinary differential equations (ODEs). Many models of biochemical systems contain 
a high number of model parameters, usually even more parameters than state vari- 
ables. A major problem in understanding biochemical systems is that most of these 
parameters are not very well known from measurements, and that they often vary sig- 
nificantly due to internal or environmental conditions of the cell. Thus analysing the 
influence of uncertain or varying parameters on stability properties is a fundamental 
issue towards understanding dynamical behaviour of biochemical systems. Moreover, 
to avoid overlooking relevant effects it is necessary to consider simultaneous changes 
in all adjustable parameters [2^[ll|. 

The requirement of looking at simultaneous changes in several parameters makes the 
application of classical continuation methods problematic, as these require to define 
a line in parameter space along which equilibrium points are tracked. A good choice 
of this line is essential to obtain meaningful results, yet this choice is often done by 
intuitive understanding of the system in the better case or iterative trials in the worse. 
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Often only a single parameter is varied at a time, but then again the choice of the 
parameter to vary is not trivial and needs to be done for example via sensitivity 
considerations. 

In this paper, we present a new method to locate points in a possibly high-dimensional 
parameter space for a change in stability properties of equilibrium points, often hinting 
to either emergence or loss of complex dynamical behaviour. The method is based on 
considering the dynamical system as a closed loop feedback system. It is then possible 
to study properties of the original system in terms of an adequately defined open loop 
system. If the open loop system is well chosen, then its dynamical behaviour is much 
simpler than that of the closed loop system. This simplification makes it possible to 
come to conclusions that could not be obtained from the closed loop system alone. In 
particular, we show how to classify parameter values where the closed loop system can 
undergo local bifurcations of equilibrium points, based on an analysis of the open loop 
system. The obtained conditions are used to develop a numerical method for searching 
parameter values that lead to a change in stability properties. In theory, this can be 
done for parameter spaces of arbitrary dimension, as neither the conditions nor the 
algorithm we use depend on the dimension of the parameter space. We consider only 
codimension one bifurcations, as they are the case that is generically encountered in 
non-linear systems. 

We make use of the fact that for stability considerations, it is sufficient to look at a 
linear approximation of the system close to the equilibrium. The linearised system is 
transformed to the frequency domain for our analysis. The use of frequency domain 
methods for bifurcation analysis has already been introduced by (author?) 1, cited 
from [15]], and relevant results have also been presented by Moiola and coworkers over 
the last decade fie','!?']. 

Several authors have also studied the problem of finding bifurcations in systems with 
many parameters using geometric tools. Based on a description of vectors normal to 
a bifurcation manifold [18| . a method to search for locally closest bifurcations from 
a given reference point was developed by |5(]. These approaches can be seen as com- 
plementary to our results. A recent application of the geometric concept to biological 
systems has been discussed by [l3]. 




Our paper is structured as follows. In Section [21 we introduce the loop breaking 
approach and provide the general tools which are necessary for our method. The 
main results are presented in Section [3l a frequency domain theorem on topological 
equivalence, an existence theorem for critical parameters and a numerical algorithm 
to search for parameters yielding a change in dynamical behaviour. Moreover, we 
shortly discuss the benefits of our approach compared to a straightforward extension 
of classical tools. As an application example, the method is used in Section |4] to 
search for possible limit cycle oscillations in an ODE model of a biochemical signal 
transduction system. 
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2 The loop breaking concept 



2.1 Problem setup 

Consider a parameter-dependent nonlinear differential equation given by 

x^F{x,p), (f) 

with a; e R", p e T' C M" and : M" xP ^ R" a smootli vector field. 

The system ([IJ is studied locally at an equilibrium point. In what follows we fre- 
quently denote 

C = ix,p)eW' xV, (2) 

where x is an equilibrium point and p a corresponding parameter. We call ^ an 
equilibrium-parameter pair of the system ([T]) in the sense that x is an equilibrium 
for the parameter p. Let be a smooth connected m-dimensional manifold of 
equilibrium-parameter pairs in R" x V, i.e. 

e 7W : = 0. (3) 

In the simplest case, there is a unique equilibrium point for each p G V, and one could 
use a function x(j)) to characterise the manifold of equilibrium-parameter pairs more 
easily. However, the approach taken here is more general and also allows to consider 
e.g. saddle-node bifurcations, where existence of a unique equilibrium for each p € V 
is not given. 



2.2 Loop breaking and closed loop eigenvalues 

Mathematically, the system ([1]) is said to contain a feedback loop if the influence graph 
of its Jacobian ^ contains a nontrivial loop [3] . Let us now assume that ([T]) contains 
a feedback loop. This assumption is not restrictive, because without a feedback loop, 
the analytical expressions for the eigenvalues in terms of parameters and the state 
variables can be taken directly from the diagonal of the (possibly permuted) Jacobian 
In this case, it is usually easy to find parameter values for a change in stability 
properties of the equilibrium points. 

In the feedback loop approach, an input-output system which corresponds to the 
original system is obtained by breaking the feedback loop. As seen from the following 
definition, the original system can be recovered by closing the feedback loop again. 

Definition 1. A loop breaking for the system ^ is a pair (/, h), where f : R" x R x 
V R" is a smooth vector field and h : R" ^ R is a smooth function, such that 

F{x,p) = f{x,h{x),p). (4) 

The corresponding open loop system is then given by the equation 

x^f{x,u,p) 
y = 
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and the closed loop system ([T]) is recovered by letting u = y. Note that there is a 
direct relation between equilibrium points in the closed and the open loop system: for 
an equilibrium-parameter pair {x,p) of the closed loop system ([T]), setting the input 
u — h{x) in the open loop system ([5]) leads to {x,p) being an equilibrium-parameter 
pair of (III). We denote u — h{x). 

To deal with the question whether different equilibrium-parameter pairs in fA can 
have different stability properties, it is reasonable to consider the linear approximation 
of the system ([T]) close to some pair ^ S Al. Only the pairs where the Jacobian 
has eigenvalues on the imaginary axis are candidate points for local bifurcations. Any 
such pair ^ is called a critical pointy and is denoted as ^c- 

The linear approximation for the open loop system ([5|) in the neighbourhood of the 
equilibrium-parameter pair ^ e is given by 

where z = x-x, T]^y~u, fi^u-u, A{£,) = ^{x,u,p), B{£,) = ^{x,u,p), 

The linear approximation of the closed loop system ([1]) can then be easily charac- 
terised as follows. 

Proposition 1. The linear approximation of the system ([T]) close to ^ G is given 
by 

z = [AiO + BiOCiO) z = A^iiO^. (7) 

Proof. This follows directly from the loop breaking definition Q and the chain rule. 

□ 

The linearised open loop system ([6]) can also be described using its transfer function, 
which is defined as 



det 



G{^,s)=C{0 {si - AiOr' BiO = 



sI-A{0 ~B{0\ 

CiO ; (8) 



det(s/- 



with the complex variable s G C. 

The following lemma is a tool to characterise eigenvalues of the closed loop system 
^ by analysing the open loop system ([5]). 

Lemma 1. sq E C is an eigenvalue of Aci{^), if and only if one of the following 
conditions holds: 

(i). So is not an eigenvalue of A{^) and G(^, sq) = 1; 
(a). So is an eigenvalue of A{£_) and det (^^^^ff^^^^ ^^'^^1 ~ ^- 



\ CiO 

The proof is provided in the appendix. In the following, Lemma[T]is used with sq on 
the imaginary axis, to characterise critical points with the condition G{^c, so) = 1- 
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2.3 Critical frequencies and imaginary closed loop eigenvalues 

In this section, the transfer function G is represented as a complex rational function 
with real coefficients, i.e. 

where k{^) e R and q{^, s), r(^, s) are polynomials in s with real scalar functions of ^ 
as coefficients. 

Moreover, we make the following technical assumption. 

(Al) The transfer function G(^, •) does not have poles or zeros on the imaginary axis 
for any & M, i.e. 

e MVuj eR: k{£_)q{£_, juo) ^ and r(^, juo) ^ 0. (10) 

In addition, the degrees of q{S,, s) and r(^, s) in s are constant with respect to 

Starting from the premise that we are interested in stability changes produced by 
changing the characteristics of the feedback loop that was broken in @ , this assump- 
tion is usually satisfied. 

The notion of a critical frequency which is introduced in the next definition will be 
useful to compute possible eigenvalues of the closed loop system ^ on the imaginary 
axis. 

Definition 2. lOc G R is said to be a critical frequency for the transfer function 
G{^,s), tf 

G{tjLUc)eR. (11) 

Obviously, different values of f will result in different critical frequencies. For a 
specific ^, all critical frequencies are given by the solutions of the equation 

lm{qi^,juM^,-J^c))=0, (12) 

which is a scalar polynomial equation in luc, with coefficients that are real scalar 
functions of ^. 

We define the set of all critical frequencies for a specific ^ as 

n.iO - {c. e R I Im(q(e, JLo)r{^, -^ju)) = 0} . (13) 

Because only the imaginary part is considered, the polynomial in (I12p is odd. The 
following properties of the set ric(C) can then be shown easily. 

Proposition 2. For any ^ G M, the set ric(0 satisfies the conditions 

(i) e 

(a) LUc ilc(0 implies that —lUc G i7c(C)> 
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(iii) either f2c(C) = K or f2c(C) has finitely many elements. 

Note that — whenever Ini(g(^, jw)r(^, —joj)) is the zero polynomial, which 

in turn is the case whenever only the even powers of s in the polynomials q{S^^ s) and 
r(^, s) have nonzero coefficients. However, this typically contradicts assumption (Al), 
so we will not consider this case specifically. 

The relevance of critical frequencies for existence of eigenvalues on the imaginary 
axis is shown by the following result. 

Proposition 3. Assume that (Al) is satisfied. If jujc with cjc G R *s an eigenvalue of 
Acii^), then Uc E ric(0- 

Proof. By assumption (Al), jtUc is not an eigenvalue of A{^). By LemmalU we have 
G{^,juJc) = 1 and thus uic & ^dS,)- ^ 

The concept of critical frequencies can be understood intuitively when considering 
the Nyquist curve of the transfer function G(^, ju) . A critical frequency is any value lOc 
at which the Nyquist curve crosses the real axis. This is obviously a necessary condition 
for having G(^, jwc) = 1, which corresponds to the existence of an eigenvalue on the 
imaginary axis as shown in Lemma [1] Our concept is thus closely related to the idea 
of the gain margin for robustness analysis of linear control systems [24] . 

Since a variation of the equilibrium-parameter pair f influences the polynomial 
equation p^ . the set of critical frequencies ilc(C) ™E^y change significantly with In 
particular, the number of elements in r2c(0 in general needs not to be constant with 
respect to ^, which complicates the analysis. However, one can show that there is a 
minimal number of critical frequencies of the transfer function G{^, s), which depends 
on the number of open loop poles and zeros and whether they are located in the right 
or the left half-plane. To this end, define the number 



where (P-) is the number of poles of G(^, ■) and (z-) is the number of zeros of 
G(^, •) in the right (left) half complex plane. Under assumption (Al), a is constant 
with respect to ^ G A4. The number of elements in the set of critical frequencies can 
now be characterised by a. 

Proposition 4. Let a be defined by (|14p and assume that (Al) is satisfied. Then, 
for any ^ G A4, i^dS,) has at least a distinct elements, if a is odd, and at least a ~ 1 
distinct elements, if a is even. 

The proof is presented in the appendix. 

The above result is used to formulate the property of minimality for the set of critical 
frequencies. 

Definition 3. Under the assumptions of Proposition^ the set of critical frequencies 
ric(0 is called minimal, if it contains exactly a elements, where 



a ^ \p+ - p- + ~ z+\, 



(14) 




a, if a is odd 
a — 1, if a is even. 



(15) 
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This definition is applied in the second technical assumption we are going to make 
use of. 

(A2) The set of critical frequencies minimal for any ^ e A^. 

If (A2) holds, we can label the roots of the polynomial equation in a consistent 
way, writing 

ndO^HiO,^UO,---,^HO}, (16) 

where the luI are continuous functions of the equilibrium-parameter pair ^ and can be 
identified with different solution branches of the polynomial equation p2p . 

Given a transfer function G(^, ■) and a corresponding set of critical frequencies r2c(^). 
Proposition m can be used to easily check the minimality of ric(0- Graphically, a 
sufficient condition for minimality of ric(0 that the Nyquist curve G(^, jw) encircles 
the origin monotonically as to varies from —00 to 00. 

3 Main results 

3.1 Topological equivalence of equilibria 

Changes in stability properties of equilibrium points are most easily studied using the 
concept of topological equivalence. Here, we use a definition for hyperbolic equilibrium 
points only (see [12| for more details). 

Definition 4. Let ^1,^2 £ -M be two hyperbolic equilibrium-parameter pairs of the 
system ([T]). ^1 and ^2 are said to be topologically equivalent, if the Jacobians ^(Ci) 
and ^(^2) have the same number of eigenvalues in the left and right half-plane. 

It is a well known result from dynamical systems theory that the topological equiva- 
lence of all equilibria in two systems is a necessary condition for topological equivalence 
of the flows. Let us consider two variants of the system ([1]), one with parameter values 
pi and the other with parameter values p2 ■ In the simple case when there is only one 
equilibrium point in each variant of the system, corresponding to the pairs ^1 and ^2, 
topological equivalence of ^1 and ^2 is a necessary condition for topological equivalence 
of the flows. In applications, we are often interested in flnding parameter values p2 
such that the system ([1} changes its dynamical behaviour when varying parameters 
from initial values pi to p2. For this problem, it is sufficient to find pairs ^1 and ^2 
which are not topologically equivalent. Due to the continuous dependence of eigenval- 
ues on parameters, this can only happen when the Jacobian ^{£,c) has eigenvalues on 
the imaginary axis for some critical point ^ At this point, we can make use of 
the methodology developed in the previous section. 

To this end, consider the set of critical frequencies ^c{£,) for a specific value of the 
equilibrium-parameter pair ^. Define the number f3{^) to be the number of elements 
LOc in ric(C) such that G{£_,jujc) > 1, i.e. 

PiO = card {u;, £ n,{0 \ G{^, juo,) > 1} , (17) 
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where cardiS denotes the number of elements in the set S. 

Geometrically, if r2c(0 minimal, /3(^) gives the winding number of the graph of 
G{^,ju!) around the point 1 in the complex plane (see Lemma [3] in the Appendix). 
The argument principle can then be used to characterise topologically equivalent 
equilibrium-parameter pairs of the system ^ via the number /3(^). We first give 
some intermediate results as Lemmas before presenting the main theorem. The Lem- 
mas are proven in the appendix. 

Lemma 2. // the set of critical frequencies minimal, then in the ordered 

sequence of critical frequencies a;^(i^) < a;^(^) < • • • < cj"(^), we have 

where i ~ 2, . . . , a. 

Lemma 3. Under the assumptions of Theorem [21 the winding number of the image 
of the Nyquist curve T under the transfer function G{S,i, ■), i — 1, 2, around the point 
1 is given by 

|ii;n(G(6,r),l)| 
Lemma 4. Under the assumptions of Theorem]^ we have 

\wn{G{^,,r)A)-wn{G{^2,r),l)\ = 1/3(6) - • 

Theorem 1. Assume that (Al) is satisfied. Lef^,^ e fee two hyperbolic equilibrium- 
parameter pairs of ([T]) such that Q.c{£,i) and ^c{S,2) are minimal. Then ^ and ^ are 
topologically equivalent, if and only if 

mi) = 

The proof is given in the appendix. 

3.2 Existence of marginally stable equilibria 

Let us now turn to the problem of how to find parameter values for which a change 
in stability properties of equilibria can happen. This is equivalent to searching for 
critical points ^ at which the Jacobian ^(Cc) has an eigenvalue on the imaginary 
axis. This typically means that ^ is part of a submanifold oi Ad that separates regions 
of topological equivalence, and in view of Theorem [1] there are typically equilibrium- 
parameter pairs 6 and ^ close to ^ such that P{£,i) ^ P{^2)- Equivalently, for a 
specific critical frequency w*, the transfer function value G(6 ja;*(^)) has to cross the 
value 1 when ^ varies continuously along a path from ^ to ^ . These observations are 
formalised in the following theorem. 

Theorem 2. Assume that (Al) and (A2) are satisfied. There exists a critical point 
£,c £ M such that ^{£.c) has an eigenvalue on the imaginary axis, if and only if there 
exist E A4 such that, for some i G {1, 2, . . . , a}, 

G(6,jc^:(6))<i<G(6,jc.^(6)), (18) 

where loK^) G flc{^). In that case, ±ju;^(6) is o,n eigenvalue of ^{^c)- 
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Proof. By Lemma [U our assumptions assure that a point is critical if and only if 

Necessity. Under the condition G{£^c,j^c) = li take ^1=^2= S,c and (fT5|) follows 
trivially. 

Sufficiency. Let ^1 and ^2 be such that (fT8|) holds. Connectivity of M implies that 
there is a path from ^1 to ^2 in A^. Continuity of the critical frequency oj].{£) and 
the transfer function coefhcients result in continuity of G'(f , jcL'^(^)) with respect to ^. 
This implies existence of such that G(^c, iw^(^c)) = 1 along any path from ^1 to 

6- □ 

The proof shows that the critical point is usually far from unique. It may be 
unique if Ai is of dimension one, i.e. there is only one free parameter to vary. In 
general, one will expect that there is a submanifold of critical points in Ai separating 
regions which represent different topological equivalence classes, where ^1 is an element 
of one such class, and ^2 is an element of the other class. On this submanifold, the 
bifurcations that can be encountered generically are codimension one bifurcations. 

Using Theorem [5J it is easily possible to distinguish dynamical bifurcations from 
static bifurcations. In fact, if i is chosen such that the critical frequency = is 

considered, Aci{^c) has a zero eigenvalue, which generically corresponds to a saddle- 
node bifurcation. If a critical frequency ^^(0 7^ is considered, Acz(^c) has conjugated 
imaginary eigenvalues, and one will generically get a Hopf bifurcation. 



3.3 An algorithm for a numerical parameter search 

In this section, we discuss an algorithm to search for parameter values that will lead 
to a change in stability properties of an equilibrium point. We assume that a starting 
parameter pi and a corresponding equilibrium xi are known, which we combine in the 
pair — (xi^pi) G M. It is reasonable to assume that the pair ^1 is not critical, 
otherwise it is usually straightforward to find parameter values yielding equilibrium 
points with different stability properties. Moreover, the manifold M. is assumed to 
be defined by a nonlinear equation of the form ip{£,) = 0. Often, one can directly use 
if = F, but sometimes a modification is useful to exclude some solutions if the equation 
F(x,p) = is known to have multiple solutions. The aim of the algorithm is to find 
an equilibrium-parameter pair ^2 ^ A4 such that ^2 is not topologically equivalent to 
^1 . The main theoretical basis of the algorithm is the result of Theorem [21 Thus it is 
also possible to search specifically for either static or dynamic bifurcations on a path 
from ^1 to ^2 by choosing an appropriate critical frequency. 

In order to put the problem in the framework developed in this paper, a loopbreaking 
for the system ([ij has to be defined. Then, by looking at the resulting transfer function 
G{S,i,s), possible changes in stability properties can be determined. In particular one 
has to decide whether to search for a static or for a dynamic bifurcation. This leads 
to the choice of a critical frequency which is to be considered in the algorithm. 

Denote the transfer function value for the critical frequency at a point ^ as 7(^). At 
the starting point fi, this value can be computed as 

7(ei) = G(a,j-c.^(a)). (19) 
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Note that an analytical expression of the function 7 can be derived directly, maybe with 
the support of computer algebra for more complex systems. This derivation requires 
only basic algebraic manipulations, differentiation and matrix inversion, which can all 
be done symbolically for typical system classes. In particular, it is not required to 
have an analytical solution of the equation F(x,p) = to construct 7. 
Now two cases have to be distinguished. 

1. If 7(^1) < 1, the algorithm searches a pair ^2 G such that 7(^2) > 1. 

2. If 7(^1) > 1, the algorithm searches ^2 £ such that 7(^2) < 1- 

The algorithm we are using is best described by the term gradient- directed continua- 
tion method. Continuation methods 2l| are popular in numerical bifurcation analysis, 
where they are used to trace the equilibrium curve in the combined state-parameter 
space. In our algorithm, continuation is used to stay on the manifold M. However, 
a continuation method alone is not sufficient, as M is m-dimensional with typically 
m > 1. Thus, the continuation is complemented with a gradient ascent or descent 
approach to achieve the desired value for 7(^2). 

In detail, the algorithm works as follows. We are discussing case 1 only, small 
extensions are required for dealing with both cases. 

1. Initialisation. Set =^1. Choose numerical parameters: A7 for the minimal 
required change in 7(^) per iteration, (5^*^^ as the initial step size and Smin (Smax) 
as minimal (maximal) step size. 

2. Prediction step. This step assures the desired increase in 7(^). 

a) Compute the gradient V7 

b) Compute the subspace which is tangent to M. in the point 




(20) 



c) Project V7(^(^)) onT^^M: 




(21) 



d) Set the predicted point 



(22) 



Step size control is used in the sense that (5^*' is varied to assure that the 
condition 




is satisfied, while keeping (5, 



< S'^'^ < S, 



max ■ 
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3. Correction step. Generally, ^prt^^ ^ -M, so a correction step is required to 
achieve ^(*+^) e Al. To this end, the Gauss-Newton method is used to solve the 
nonlinear equation 

^(C(^+i))=0 

for where £,pri'^^ is used as starting point for the Gauss-Newton algorithm. 

If the Gauss-Newton algorithm converges, the algorithm takes the solution as 
value for ^('+^) and proceeds to the next step. Otherwise, the algorithm reduces 
the step size (5^*' and goes back to 2d). 

4. Finishing criterion. Compute 7 (^^'^^''). If 7 (C*'*^^'') > li finish successfully, 
otherwise iterate to step 2. 

If the algorithm finishes successfully, it does so in a finite number of steps with a 
previously known upper bound due to step size control via inequality (|23|) . 

However, in the same way as classical continuation methods, the algorithm may 
fail if the Gauss-Newton algorithm in step 3 does not converge, and the step size J'-'-' 
may not be reduced further due to the constraint Smm 5: '5^*-' at the same time. This 
problem may appear if the system is numerically ill-conditioned, but can typically 
be avoided by choosing a smaller value for either A7 or for 5min, with the drawback 
of increased computational effort. Also, it can in general not be excluded that the 
function 7(^) has local extrema, which may pose problems to the algorithm. Such 
problem may be detected numerically from the vector u^*) taking very small values. 
However, in several applications we have not encountered this problem so far. 

The algorithm as described above does not consider constraints on the parameters 
p. Such constraints can be included by slight modifications in steps 2 and 3. If the 
border of the set V is approached during the iteration, the modified algorithm projects 
the gradient V7 {S,^-^^) on the intersection of the tangent to Ai and the tangent to the 
border of V. With additional step size control, a constraint violation is then avoided. 

3.4 Discussion of the feedback loop approach 

Two key steps in the approach we have taken are the transformation of the problem 
to the frequency domain and the consideration of the critical frequencies. These steps 
require an elaborate setup and therefore need to be well justified. 

Approaching the given problem in the time domain would typically require to deal 
with the eigenvalues of the Jacobian ^ on the considered manifold of equilibrium 
points. In particular, it would require to consider how the eigenvalues change if the 
parameters change. Going to the frequency domain will typically reduce the num- 
ber of variables that are to be tracked with changes, because there are typically less 
critical frequencies than eigenvalues. For a minimum phase system, the number of crit- 
ical frequencies is not more than the relative degree of the transfer function G(^, s). 
Moreover, the position of eigenvalues has to be tracked in the two dimensions of the 
complex plane, whereas the transfer function values at critical frequencies are always 
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real numbers. In particular, it would be difficult to estimate from the eigenvalues of 
the system for some starting parameters, which pair of eigenvalues should be pushed to 
the imaginary axis in order to obtain a Hopf bifurcation. Using the frequency domain 
approach, the critical frequency for which the transfer function value should be pushed 
towards 1 can typically be determined easily. 

In classical bifurcation analysis, so called bifurcation test functions are used to check 
whether a bifurcation may occur when going from one parameter value to another one 
[l3 |. The test function ^E* is defined such that "^{^c) = if the bifurcation that is 
tested for occurs at ^c- Bifurcations are detected by the test function ^'(^) changing 
sign when going from one point to the other, i.e. if 5'(^i)^'(C2) < 0, then a bifurcation 
occurs between and ^2- For bifurcations of codimension one, suitable test functions 
are known and are routinely used in numerical continuation algorithms. Note that in 
the frequency domain approach, the expression G'(^, jWc(C)) — 1 is a test function for 
a generic saddle-node bifurcation, if we consider ujc — 0, and it is a test function for a 
generic Hopf bifurcation when considering uJc ^ 0. Computing classical test functions 
for a given point ^ requires a similar or slightly less computational effort as computing 
the transfer function values at the critical frequency. So we need to justify why we do 
not use classical test functions for bifurcation search in a high-dimensional parameter 
space. 

Because classical continuation methods cannot be used in a high-dimensional pa- 
rameter space, one has to look for different approaches. A naive approach to find 
parameters for a bifurcation would be to solve directly the equations 



However, in most cases this will be numerically infeasible with classical test functions, 
even if the combined parameter state space is of very low dimension. A more sophis- 
ticated approach could use basically the same algorithm that we have presented in 
Section 13.31 the gradient-directed continuation method, and just use the gradient of 
a classical bifurcation test function instead of the gradient of the transfer function 
juJc{£,))- We have also implemented this approach for several examples, but run 
into numerical problems for any system of medium complexity. In particular, the anal- 
ysis presented in the next section did not work with a classical bifurcation test function 
for a Hopf bifurcation due to numerical problems. These problems seems to be related 
to our observation that the value of the classical bifurcation test function seems to 
be numerically much less well behaved with respect to parameter variations than the 
transfer function value at critical frequencies. This observation is also supported by 
the fact that general optimisation methods can indeed find suitable parameter values 
for a bifurcation by directly solving (PS)) . if the expression G{£_, jcudO) — 1 is used as 
a test function \2M . 
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Figure 1: Illustration of the MAPK cascade model with reaction numbers. 



4 Application to a biological signalling system 

In this section, we apply the theoretical results and the numerical algorithm described 
in the previous section to a system for biochemical signal transduction. A central 
element of the signal transduction in eukaryotic cells is the mitogen activated protein 
kinase (MAPK) cascade. It appears in several signalling pathways and is related to 
cell differentiation, proliferation, and response to external stress 30l 20| . Several ODE 
models for this system have been proposed during the last decade [l9| . 

The MAPK cascade consists of three layers of kinase proteins, where each kinase 
activates the next layer, and the last layer corresponds to the output of the cascade. 
We consider the MAPK cascade as it appears in the EGF (epidermal growth factor) 
receptor pathway 0, [1^ . There, the three kinases in the order how activation proceeds 
are Raf, MEK (MAPK/ERK kinase), and ERK (extracellular-signal-regulated kinase) . 
Active ERK phosphorylates and inhibits SOS (son of sevenless homologue), which is 
required in the activation of Raf [1] . This constitutes a negative feedback loop in the 
system. The model we use here is a slight simplification of a model suggested by 
and it is also a subsystem of the EGF pathway as modelled by [3] . A cartoon of the 
biochemical reactions incorporated in the model is shown in Fig. [TJ 

In the equations, the concentrations of phosphorylatcd kinases are denoted as xn = 
[Raf*], X21 = [MEK-P], X22 = [MEK-PP], X31 = [ERK-P] and 2:32 = [ERK-PP]. The 
concentrations of unphosphorylated inactive kinases Raf, MEK and ERK need not be 
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Table 1: Reaction rates in the MAPK cascade model 



Reaction 


Rate 


Vl 


T/ 2^1t — ^11 


V2 




V3 


k3Xll{x2t - X21 - X22) 


Vi 


^4X11X21 


V5 




V6 




V7 


hx22{xzt - 2^31 - 3=32) 


VS 


k%X22X3\ 


V9 


\r X32 

K„,g+X32 


VW 


^'^^K^^o+xs^ 



included as state variables, as they can be computed via the conservation laws 

[Raf] + xii = xit 
[MEK] + X21 + X22 = X2t 
[ERK] + + X32 = X3t, 

where xu, X2t and 2:34 are parameters for the total concentrations of kinases, which 
are constant. Table [1] shows the mathematical expressions for the reaction rates, with 
numbers corresponding to the labels in Fig. [1] Nominal parameter values for the 
simplified model have been adopted from [10| , and are shown in Table [2] as pi . 

Using the reaction rates from Table [H the model can be written as a system of five 
ODEs with 20 parameters: 



Xu 


= Vl 


- V2 






X21 


= V3 


+ V5 


- V4 


- Ve 


X22 


= V4 


- V5 






X31 


= V7 


+ Vg 


- Vs 


- Vio 


X32 


= 1^8 


- Vg. 







The only difference to Kholodenko's model is in the phosphorylation reactions 3, 
4, 7 and 8. The original model uses Michaelis-Menten kinetics for these reactions, 
whereas our simplified model uses mass action kinetics. It can be argued that with the 
concentrations of all kinases being on the same order of magnitude, the assumptions 
for using Michaelis-Menten kinetics in reactions 3, 4, 7 and 8 are not valid anyway, 
and one could aim to achieve a similar dynamical behaviour with the simpler model 
structure where mass action kinetics are used. 
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Kholodenko has shown in simulations that the system can show hmit cycle oscilla- 
tions for some parameter values. Due to the simplifications in four reaction rates, the 
model docs not oscillate for nominal parameter values pi. Instead, the model has 
a stable equilibrium xi for these parameter values. Solutions of the model converge to 
the steady state within 20 seconds, as depicted in Fig. [3l 

The question we deal with is whether parameters can be changed such that also the 
simplified model shows sustained oscillations. To answer this question, we apply the 
algorithm to search for critical parameter values which is described in the previous 
section. 

The first step in our analysis is to choose a suitable loop breaking. For the given 
system, an intuitive approach is to break the loop at the feedback inhibition of reaction 
vi by ERK-PP. Thus we choose h{x) = X32 and replace X32 by the input u in the 
reaction rate wi, thus obtaining the dynamics of the open loop system f{x, u,p). 

A linearisation of the open loop system around the equilibrium point and a Laplace 
transformation give the transfer function G(^, s), whose graph is shown in Fig. [2] The 
problem is now to find parameters p2 with a corresponding unstable equilibrium point 
X2- This can be done using the numerical algorithm presented in Section [3.31 

The set of critical frequencies is minimal with a = 3, which can be seen from Fig. [2] 
by the observation that the graph of G(^i, jw) encircles the origin monotonically and 
crosses the real axis three times. The origin of the complex plane is not counted as 
crossing, as we have cu = ±cxd there. The only positive critical frequency is w^(^i) = 

0. 017 s""'^, and we will consider only lo^ in the search for destabilising parameters, 
because our goal is to find a Hopf bifurcation. The corresponding transfer function 
value is G{£,i,jujl{^i)) — 7(^1) = 0.12, corresponding to the equilibrium xi being 
stable in the closed loop system. 

The goal for the parameter search algorithm is to find parameters such that 7(^2) > 

1. Then the corresponding equilibrium point X2 will not be topologically equivalent 
to the nominal equilibrium xi and we can expect a Hopf bifurcation when varying 
parameters from the nominal value pi to the new value p2- To ensure that the algorithm 
does not stop at the bifurcation point, but continues to vary parameters until the 
oscillations have reached a considerable amplitude, we try to achieve 7(^2) > 1-5 in 
the implementation used here. 

In the application of the algorithm to this problem, we set the minimal change in the 
transfer function value per iteration A7 = 10~^ and the initial step size S^*^^ = 10"'*. 
The Gauss-Newton algorithm in the correction step was constrained to 20 iterations, 
but in the step size control the step size ^^^'^ was already decreased if the Gauss-Newton 
algorithm required five or more iterations for convergence. With these settings the 
algorithm finishes successfully after 276 iterations, yielding the parameters p2 and an 
equilibrium X2 with the transfer function value G{^2, j^ci^'i)) = 1-52 and the critical 
frequency ^^(^2) = 0.0068 s""'^, where ^2 = (^2,^2)- The parameter values in p2 
are shown in Table O Note that although in principle all parameters could have been 
changed when going from pi to p2 , the algorithm varies only 9 out of the 20 parameters 
by an amount of more than 20 %. Since the algorithm uses the gradient of the transfer 
function value at the critical frequency, we can presume that the parameters that have 
been varied by a larger amount have higher influence on existence of oscillations than 
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Figure 2: Nyquist plots of open-loop MAPK model for parameters pi (grey line) and 
P2 (black line). 

the other parameters. 

The graph of G{£,2,j^) is shown in Figured) For the new parameters p2, the graph 
now encircles the point 1. By Theorem [TJ we see that the equilibria xi and X2 are not 
topologically equivalent. Indeed, X2 is unstable and the system converges to a limit 
cycle for parameters p2- The time course of these oscillations is plotted in Figured 

In conclusion, our method is able to compute parameters which render the corre- 
sponding equilibrium unstable and thus leads to the emergence of sustained oscilla- 
tions. About half of the parameters are varied by a non-negligible amount, but all 
variations are within the physiological range, the highest variation being a factor of 
about 15 in the Km-yalue of one reaction. It is also worth mentioning that the con- 
centration values in the equilibrium did not change significantly, although this should 
not be of physiological relevance for the unstable equilibrium. 

5 Conclusions 

The loop breaking concept is introduced as a theoretical tool to analyse complex 
behaviour in ODE systems as frequently encountered in mathematical biology. Based 
on this tool, we present results on topological equivalence of equilibria in systems with 
high-dimensional parameter spaces and on the existence of critical parameters, for 
which stability properties of equilibria may change. In addition, an algorithm is given 
to systematically search for critical parameters. Using an ODE model for a MAPK 
cascade, we show that the algorithm can be used to efficiently search for parameter 
values leading to limit cycle oscillations in the system. 

Non-uniqueness of critical parameter values is a problem that is inherent to this kind 
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Table 2: Reference parameters pi and parameters for instability P2 in the MAPK 
cade model. 



Par am. 


Pi 


P2 


Unit 


rel. change 




2.5 


2.5 


nM/s 


1.00 




9 


18.9 


nM 


2.09 


Kml 


10 


8.1 


nM 


1.23-1 


V2 


0.25 


0.17 


nM/s 


1.43-^ 


Km2 


8 


0.54 


nM 


14.8-1 


k3 


0.001 


4.3 • 10-"* 


l/(s nM) 


2.34-1 


k4 


0.001 


5.7-10^* 


l/(s nM) 


1.76-1 




0.75 


0.74 


nM/s 


1.01-1 




15 


7.6 


nM 


1.98-1 




0.75 


0.77 


nM/s 


1.02 




15 


13.9 


nM 


1.08-1 


k7 


0.001 


5.2 • 10""* 


l/(s nM) 


1.93-1 


ks 


0.001 


7.9 • 10-4 


l/(s nM) 


1.27-1 


V9 


0.5 


0.49 


nM/s 


1.02-1 


Krn9 


15 


15.1 


nM 


1.01 


Vio 


0.5 


0.51 


nM/s 


1.01 




15 


15.4 


nM 


1.03 


Xlt 


100 


100.2 


nM 


1.00 


^21. 


300 


300.2 


nM 


1.00 


•*:!/, 


300 


304.4 


uM 


1.01 
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Figure 3: Convergence to steady state for parameters pi (grey line) and sustained 
oscillations for parameters p2 (black line). The oscillations coexist with the 
unstable equilibrium X2 (dotted line). 

of analysis. If the dimension of the parameter space is higher than the codimension of 
the bifurcation, then there will be a submanifold of bifurcation points in the parameter 
space. Our algorithm computes one of these points. Starting from a critical point thus 
found, one can then use continuation methods to further explore the structure of the 
set of critical points. 

Another possibility for further studies would be to search for a bifurcation which 
is locally closest to some reference parameter values. A method for this has been 
presented by Q . The method requires a bifurcation point where the search is started, 
and we expect our algorithm to give a starting point which is better suited for the 
method discussed in [5| than a bifurcation search along a random line in parameter 
space. 

The biological example we study in Section [4] is simple in that it contains only one 
feedback loop. Of course many biological systems contain several intertwined feedback 
loops. Then, the choice of the loop breaking point needs more attention. However, 
a comparison of different loop breaking points could give hints to the role played by 
individual loops in the dynamical behaviour of the system. 
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Appendix 
Proof of Lemma [1] 

For simplicity of notation, wc drop the dependence on ^ of matrices A, B and C. By 
Schur's lemma, we have 



det(s/ -A- BC) = det 



si -A -B 
-C 1 



Let (sJ- e m("-i)x" denote the matrix {sI~A) with the i-th row deleted. Then, 
by cofactor expansion, 

si -A -B\_^ ^^^^^^ _,^„+i+,^ /^(s/- 

-C 



det y _^ 1 J " ^ ■ '^^^^''^ ^{-^T^^^'h det 

In the same way, 

det (^^_"/ ^ - det ((^^ 



and with Prop. [T] it follows that 

det(s/ - Aci) = det(s/ -A)- det (^^ ~ ^ . (27) 

So is an eigenvalue of Ad if and only if det(so/ — Ad) = 0. For condition (i), we 
have det(so-^ ^ ^) 7^ 0, and thus the equation 



det 



(sqI ~ a -B 
\ C 
det(so/ - A) 



is equivalent to Sq being an eigenvalue of Ad- The claim then follows from ([5]). 

The other case where det(so/ — ^) = is considered in condition (ii), and (|27p can 
be used directly to prove the claim. 



Proof of Proposition [4] 

Note that a is constant over M. due to assumption (Al). Consider the transfer function 
G'(^, s) for a constant ^ G A^. For ease of notation, we drop ^ in the transfer function 
in the following. 

It is well known from linear control theory that the argument of G{jui) changes by 
air when varying to from — oo to oo 4]: 

I argG(joo) — argG(— joo)| — air. 
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The symmetry G(jw) = G{—juj) implies that argG(joo) = — argG(— joo). From these 
two facts, it follows that the argument of G{juj) spans the open interval la — (—^7 ^) 
for uj G (—00, 00). 

Moreover, by Definition [2] the condition ujc G ^c{£,) is equivalent to 

a.igG{jujc) = kn, k eZ. 

If a is even, the claim follows directly, since there are a — 1 different integer values 
for k such that kir is inside the interval . This corresponds directly to having a — 1 
or more critical frequencies. 

In the other case, if a is odd, some additional reasoning is needed to prove the 
proposition. In this case one has 

2m + 1 2™ + 1 \ 

TT, TT , m G Z. 

2 '2 y ' 

Thus the borders of the interval are not at integer multiples of tt, which implies that 
in this case there are a different integer values for k such that kn is in ; corresponding 
to at least a critical frequencies. 

Proof of Theorem [T] 

The proof for Theorem [1] uses the^argumcnt principle from complex analysis, which is 
repeated here for completeness |2t^ . 

Theorem (The argument principle): Let f be a meromorphic function on the 
domain D <Z C and T a simply closed curve in D such that f does not have a zero 
or pole on T. The winding number wn(f{T), 0) of the image of T under f around the 
origin is given by 

wn{f {r),0) = Zf -pf, 

where Zf (pf) is the number of zeros (poles) of f in the interior of the curve T , counted 
according to their algebraic multiplicities. 

Note that the winding number is counted in the counter-clockwise direction. 

As typically done in linear control theory, we will generally use the imaginary axis 
for r, also called the Nyquist curve. This can be seen as a closed curve by first taking 
only the interval [—jR, jR] and the half circle with radius R in the right half plane, 
and second letting R — > 00. Thus the interior of T is the right half plane. 

We will first proof the intermediate results given in Lemmas UHll before proving 
Theorem [11 

Proof. (Lemma\^ If ric(C) is minimal, then there is exactly one Uc S f^c(C) such that 
argG(^, jciJc) — kn for each fc G Z with kn G la (where la is the interval defined in 
the proof of Proposition 2]) . This implies that 

|argG(e,K(C)) - argG(e,K"'(e))| = ^ 

and thus G{^,ju;l{0)G{tjcol-\0) < 0. □ 
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Proof. (Lemma\B^ Note that the loop breaking ^ assures that G(^, joo) = —joo) = 
0. Considering also Lemma [21 it follows that every cut of G{^i,r) to the right of the 
point 1 is preceded and followed by a cut of G{£_i,r) with the negative real axis. Thus 
each cut to the right of the point 1 corresponds to one winding of G(^i, F) around the 
point 1. Moreover, Lemma [2] assures that these windings all have the same direction 
and thus several windings cannot cancel in the total winding number. □ 

Proof. (Lemma^ From assumption (Al), the transfer functions G(^i, •) and G(^2, •) 
have the same number of zeros and poles in the left and right half plane. Thus for the 
phase differences we have 

argG(^i,joo) - argG(Ci, -joo) = arg G(C2, joo) - arg G(6, -joo). 

This implies that the winding numbers z«n(G(^i, F), 1) and z«n(G(^2, F), 1) have the 
same sign and with Lemma |3] we conclude 

\wn{G{^,,r), 1) - u;n(G(6, F), 1)| = | |u;n(G(Ci, F), 1)| - |u;n(G(e2, F), 1)|| 

= |/3(a)-/3(6)|- 

□ 

We are now ready to give the proof of Theorem [T] 

Proof. {Theorem]!^ By Definition 31 topological equivalence of and ^2 is equiva- 
lent to the condition that the matrices ^c/(Ci) ^^^d (^2) have the same number of 
eigenvalues with positive real part. 

From the proof of Lemma [1] we know that 

detjsl ~ A^iiO) 
^^^'')^^- det(s/-A(0) ■ 

Using the argument principle, it follows that 

where ncz(^) {rioiiC)) is the number of eigenvalues of Aci{£,) {A{£^)) with positive real 
part. By assumption, rioii^i) — noi(S,2) and thus ^1 and ^2 are topologically equivalent 
if and only if 

«;n(G(ei,F),l) -ii;n(G(C2,F),l). 
The claim of the theorem then follows from Lemma HI □ 
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